    # rm(list=ls()) ; gc()
    
    
    #Article:    The Politicization of Bureaucrats: Evidence from Brazil			
    #Authors:     Anderson Frey <anderson.frey@rochester.edu>	& Rogerio Santarrosa <rogeriobs2@insper.edu.br>		    
    
    # This script produces all Tables and Figures in the Online Appendix
    # Tables are produced in a way that they can be directly copied/pasted in a LaTex document
    
    # Produced and replicable with RStudio v2024.04.1 and R v4.2.1

    # It needs the following files:  data_main.csv (Tables B1,C2-C7,C9,C10,C12, E1 & Figures C1-C3 ) 
                                    #data_cadunico.csv (Tables B1,C1,C2,C5,C6,C8,C11)
                                    #data_bureau.csv (Tables C1,C5,C6,C8,C9)
                                    #data_partisans.csv (Tables C5,C6)
                                    #data_matched.csv (Table C8)
    
    # Please upload all libraries and functions below before attempting to produce any Tables or Figures

    # The seed for all bootstraps in set at the University of Rochester ZIP CODE: 14627   
 
    # Set the working directory to the folder that contains the replication files, use the function below

    setwd("")
    
    
    # Install Packages ==============================================================================================================
    
      install.packages("stringr") 
      install.packages("lubridate") 
      install.packages("data.table") 
      install.packages("dplyr") 
      install.packages("stringi")
      install.packages("lfe")
      install.packages("rdrobust")
      install.packages("fastDummies")
      install.packages("ggplot2")
      install.packages("boot")
      install.packages("gridExtra")
    
    # ==============================================================================================================================
    
    
    # Libraries ====================================================================================================================
      
      library(stringr) 
      library(lubridate) 
      library(data.table) 
      library(dplyr) 
      library(stringi)
      library(lfe)
      library(rdrobust)
      library(fastDummies)
      library(ggplot2)
      library(boot)
      library(gridExtra)
    
    # =============================================================================================================================
    
    
    # Functions ===================================================================================================================
    

      .LatexTable <- function(data,rownames,spaces,dist){
        cols <- ncol(data)*2+1
        res <- data.frame(matrix(0,nrow(data),cols))
        res[,] <- "&"
        fill <- seq(from=2,to=cols,by=2)
        res[,fill] <- data
        res <- cbind(rownames,res)
        lastone <- rep("\tabularnewline",nrow(res))
        lastone[spaces] <- paste0("\tabularnewline[",dist,"cm]")
        lastone[nrow(res)] <- paste0("\tabularnewline[",dist,"cm]")
        res[,ncol(res)] <- lastone
        names(res) <- NULL
        return(res)
      }
      
      
      .RDPlot <- function(data,cut,band,bins,degree,x.title,y.title,title,con,sz,brk,limi1,limi2) {
        
        if(degree == 1){mm <- (var ~ Score2 )}
        if(degree == 2){mm <- (var ~ Score2 + I(Score2**2) )}
        if(degree == 3){mm <- (var ~ Score2 + I(Score2**2) + I(Score2**3) )}
        
        data <- subset(data, !is.na(var))
        
        data$Score2 <- data$run - cut 
        data <- subset(data, abs(Score2) <= band)
        vec <- seq(cut, band, length.out =  (bins+1))
        vec_mid <- vec[-length(vec)] + diff(vec)/2
        data$count <- 1 ;     i <- 0
        
        
        for (iii in 0:1){
          cc <- iii 
          if(iii == 0) cc <- -1
          aux <- subset(data, t == iii)
          aux$groups <- as.factor(cut(cc*aux$Score2,breaks=vec, include.lowest = TRUE))
          R <- group_by(aux,groups) %>% summarise (var = mean(var), count = sum(count) ) #%>% complete(groups)
          R$t <- iii
          R$Score2 <- cc*vec_mid
          
          mdl <- lm(mm , data = aux) 
          pred <- predict(mdl, R, interval = "confidence", level = con, na.action = na.pass) 
          PRE <- pred[,1] ; R$pre <- PRE
          predSE <- predict(mdl, R, interval = "confidence", level =con, na.action = na.pass, se=TRUE)$se.fit
          t <- ((pred[1,1]-pred[1,2])/predSE[1])
          V <- as.matrix(vcov(mdl) )  
          
          if(degree == 1){X <- data.frame(1,R$Score2)}
          if(degree == 2){X <- data.frame(1,R$Score2,(R$Score2^2))}
          if(degree == 3){X <- data.frame(1,R$Score2,(R$Score2^2),(R$Score2^3))}
          
          X <- as.matrix(X) ; se_robust <- sqrt(diag(X %*% V %*% t(X)))
          LOW <- pred[,1] - t*se_robust ; HIGH <- pred[,1] + t*se_robust
          R$low <- LOW ;    R$high <- HIGH 
          
          if(iii == 0){R0 <- R ; plot_data <- R}else{R1 <- R ;plot_data <- rbind(plot_data,R)}
        }
        
        plot_data <- data.frame(plot_data)
        plot_data <- filter(plot_data, !is.na(count))
        plot_data$count2 <- plot_data$count/sum(plot_data$count)
        if(sz==0)  plot_data$count2 <- 1
        ymx <- max(plot_data$var)
        ymn <- min(plot_data$var)
        
        ggplot(data=plot_data, aes(x=Score2, y=var, ymax = ymx, ymin = ymn))   +   theme_bw() +
          theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
          ylim(limi1,limi2)+
          geom_point( data=R0, aes(x=Score2, y=var, size=count), shape=21,  color = "gray25", bg="gray45") + #aes(size = count2/4)
          geom_point( data=R1, aes(x=Score2, y=var,size=count), shape=21,   color = "gray55", bg="gray75") +
          geom_vline(xintercept = 0, colour="gray25", linetype = "longdash") +
          geom_line(data=R1, aes(x=Score2, y=pre),colour="gray55",  linewidth=2, alpha=.95) +
          geom_line(data=R0, aes(x=Score2, y=pre),colour="gray25",  linewidth=2, alpha=.95) +
          #ggtitle(title) + 
          #theme(plot.title = element_text(family = "A", size=12))+
          xlab(x.title) + ylab(y.title) +
          theme(axis.text = element_text(colour="black", size=12, family="A")) +
          theme(axis.title = element_text(colour="black", size=12, family="A")) +
          geom_ribbon(data=R0, aes(ymin=low,ymax=high),alpha=0.10) +
          geom_ribbon(data=R1, aes(ymin=low,ymax=high),alpha=0.10) +
          scale_x_continuous(breaks=brk) +
          theme(legend.position="none")
        
        
      }
      
      
      .DiffRDD1 <- function( cities, base, indices, b1a, b1b, v1, v2, con, byvar){
        
        library(dplyr)
        library(data.table)
        library(lfe)
        
        b <- data.table( cities[indices,] )
        bb <- data.table( inner_join(b,base,by="ibge") )
        
        bbx <- filter(bb, abs(run) < b1a )
        if(byvar==1) bbx <- filter(bbx, i==0 )
        bbx$var <- bbx %>% select(all_of(v1))
        if(con==0) reg <- felm(var ~ t*run  |0|0|0 ,  data=bbx) 
        if(con==1) reg <- felm(var ~ t*run +age+gender+public_servant+newcomer+education+past_memb1+past_memb2+past_memb3+pt_coalition+pt_mayor+pt_federal+size  |0|0|0 ,  data=bbx) 
        if(con==2) reg <- felm(var ~ t*run +age+gender+public_servant+newcomer+education+past_memb1+past_memb2+past_memb3+size  |0|0|0 ,  data=bbx) 
        a <- summary(reg)$coefficients["t",1]
        
        bbx <- filter(bb, abs(run) < b1b )
        if(byvar==1) bbx <- filter(bbx, i==1 )
        bbx$var <- bbx %>% select(all_of(v2))
        if(con==0) reg <- felm(var ~t*run |0|0|0 ,  data=bbx) 
        if(con==1) reg <- felm(var ~t*run +age+gender+public_servant+newcomer+education+past_memb1+past_memb2+past_memb3+pt_coalition+pt_mayor+pt_federal+size |0|0|0 ,  data=bbx) 
        if(con==2) reg <- felm(var ~t*run +age+gender+public_servant+newcomer+education+past_memb1+past_memb2+past_memb3+size |0|0|0 ,  data=bbx) 
        b <- summary(reg)$coefficients["t",1]
        
        return(a-b)
        
      }
      
    # =============================================================================================================================
   
      
    # Basic Setup | Run before starting Tables and Figures =========================================================================
      
      windowsFonts( A=windowsFont("Liberation Sans")) # Font for ggplot2 
      
      party_list1 <- c(15,22,40,12,65,10,20,19,36,13) # Parties in PT's federal base in 2010
      party_list2 <- c(15,22,40,12,65,10,20,19,36)    # Parties in PT's federal base in 2010 (ex-PT)
      
      controls <- c("age", "gender", "public_servant", "education", "newcomer","pt_mayor", "pt_federal", "pt_coalition", "size", # define covariates
                    "past_memb1", "past_memb2", "past_memb3")
      

      btype <- "cerrd"  # type of algorithm for bandwidth calculation
      qq <- 2           # polynomial for bandwidth calculation
      vtype <- "hc3"    # type of standard error for the RD estimation
      
    # =============================================================================================================================
  
      
    # Necessary Files | Run before starting Tables and Figures=====================================================================
      
      cadunico <- fread( "data_cadunico.csv")
      MainData <- fread( "data_main.csv")
    
      
    # =============================================================================================================================
      
  
    # Table B1 | Patronage | Old members | Old Interviewers ========================================================================
          
        w1 <- w2 <- MainData
        w1$i <- 0
        w2$i <- 1
        wxx <- rbind(w1,w2)
        
        
        newzy <- cadunico
        newzy$t <- newzy$coalition_side_past
        
        newzy$inter_tot <- 1 
        newzy$inter_rel <- ifelse(!is.na(newzy$old_party),1,0)
        newzy$i <- ifelse(newzy$date_int <= "2008-09-30"  ,1,0 )
        
        
        aa <- data.table(group_by(newzy,ibge,t,i) %>% summarise(outcome.new=sum(inter_rel),inter_tot=sum(inter_tot)))
        bb <- data.table(group_by(aa,ibge,i) %>% summarise(total.new=sum(inter_tot)))
        aa$inter_tot <- NULL
        
        wxx <- left_join(wxx,bb,by=c("ibge","i"))
        
        wxx <- left_join(wxx,aa,by=c("ibge","t","i"))
        wxx[is.na(wxx)] <- 0
        
        wxx$var <- wxx$outcome.new*100 / wxx$total.new
        wxx <- filter(wxx, !is.na(var))
        
        bands <- rdbwselect(wxx$var, wxx$run, c = 0, p = 1, q = qq,bwselect = btype , masspoints="off")$bws ; bands*.1
        b1 <- bands[1] ; b1
        b2 <- bands[3]
 
                
          res <- data.frame(matrix(0,5,3))
      
          for (j in 0:1){
            
            wx <- filter(wxx, i == j)  
             
            if(j==0) b1a <- b1
            if(j==1) b1b <- b1
            
            covars <- wx %>% select(all_of(controls))
             
            reg <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars,vce=vtype) 
            reg2 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars, level=90, vce=vtype) 
            reg3 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, vce=vtype) 
            
            if(j==0) c1 <- reg$coef[1]
            if(j==1) c2 <- reg$coef[1]
            co <- round(reg$coef[1],3)
            se <- round(reg$se[3],3)
            ci1 <- round(reg$ci[3,1],3)
            ci2 <- round(reg$ci[3,2],3)
            n <- round(reg$N_h[1],0)
            stars <- ""
            if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
            if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
            bas <- round(reg3$tau_cl[1],3)
            if(j==0) bas1 <- bas
            if(j==1) bas2 <- bas
            
            res[1,j+1] <- paste(format(co,nsmall=3),stars,sep="")
             res[2,j+1] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
            res[3,j+1] <- paste(format(bas,nsmall=3),sep="")
            
            dt <- filter(wx, abs(run) < b1 )
            res[4,j+1] <- length(unique(dt$ibge))
            res[5,j+1] <- format(round(b1,2),nsmall=2) 
          }
      
          
          cities <- wxx
          cities <- unique(cities, by="ibge")
          cities$x <- 1
          cities <- cities %>% select(ibge,x)
          
          set.seed(14627)
          results <- boot(data=cities, base=wxx, statistic=.DiffRDD1 , R=500,con=1,v1="var",v2="var",
                          b1a=b1a, b1b=b1b,  parallel = "snow", ncpus = 4, byvar=1)
          
          ci = boot.ci(results, type=c("perc","basic"), conf=.95, vce=vtype)  
          ci2 = boot.ci(results, type=c("perc","basic"), conf=.90, vce=vtype) 
          
          stars <- ""
          if(ci2$perc[4]*ci2$perc[5]>0) stars <- "+"
          if(ci$perc[4]*ci$perc[5]>0) stars <- "*"
          
          ci1 <- round(ci$perc[4],3)
          ci2 <- round(ci$perc[5],3)
          
          res[1,3] <- paste(format(round(c1-c2,3),nsmall=3),stars,sep="")
          res[2,3] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
          res[3,3] <- paste(format(round(bas1-bas2,3),nsmall=3),sep="")
          res[4:5,3] <- "-"
          res1 <- res
      
      
          rw <- c("Incumbency Effect","C.I. (95%)","Baseline","Municipalities","Optimal Bandwidth")
          a <- .LatexTable(res1,rw,spaces=c(2,5),dist="0.2")
          print(a, row.names=FALSE)
          
          rm(wxx,w1,w2,reg,res,bands,b1,b2,covars,newzy,aa,bb,reg2,reg3,res1,results,ci,cities,dt,a,wx)
          
          
          
    # =============================================================================================================================
        
     
    # Table C1 | Summary Statistics: CadUnico Interviewers and Other Bureaucrats ==================================================

        bureau <- fread("data_bureau.csv")
        summa1 <- filter(cadunico,  !is.na(ocup) & salary>0 & education>=0   )
        summa2 <- filter(bureau, !is.na(ocup) & salary>0 & education>=0   )
     
        summa1$ocup2 <- as.numeric(substring(summa1$ocup,1,3))
        summa1$ocup2[summa1$ocup2 %in% c(231,232,233,234)] <- 231
        summa1$ocup2[summa1$ocup2 %in% 411:415] <- 411
        
        summa2$ocup2 <- as.numeric(substring(summa2$ocup,1,3))
        summa2$ocup2[summa2$ocup2 %in% c(231,232,233,234)] <- 231
        summa2$ocup2[summa2$ocup2 %in% 411:415] <- 411
        
        
        agg <- group_by(summa1, ocup2) %>% tally()
        agg <- agg %>% arrange(desc(n))
        ocs <- agg$ocup2[1:6]
        
        
        for (i in 1:7){
          if(i < 7) job <- ifelse(summa1$ocup2 == ocs[i],1,0)
          if(i < 7) job2 <- ifelse(summa2$ocup2 == ocs[i],1,0)
          if(i == 7) job <- ifelse(!summa1$ocup2 %in% ocs,1,0)
          if(i == 7) job2 <- ifelse(!summa2$ocup2 %in% ocs,1,0)
          summa1 <- cbind(summa1,job)
          summa2 <- cbind(summa2,job2)
          names(summa1)[ncol(summa1)] <- paste0("job",i)
          names(summa2)[ncol(summa2)] <- paste0("job",i)
          
        }
        
        
        summa1$college <- ifelse(summa1$education > 8 ,1,0)
        summa2$college <- ifelse(summa2$education > 8 ,1,0)
        
        summa1$age <- 2008-summa1$year_born 
        summa2$age <- 2008-summa2$year_born 
        
        
        for(i in 1:2){
          
          if(i==1) xc <- summa1
          if(i==2) xc <- summa2
          
          xc$par <- ifelse(!is.na(xc$old_party),1,0)
          xc$par[!is.na(xc$new_party) & !is.na(xc$old_party)] <- 2
          xc$par[!is.na(xc$new_party) & is.na(xc$old_party)] <- 3
          xc$count <- 1
          
         
          
          if(i==1) agg <- group_by(xc, par) %>% summarise( count=sum(count)/1000, 
                                                           wage=round(mean(salary),0), age=mean(age),
                                                           college=mean(college), female=mean(gender),newbie=mean(year_hired),
                                                           job1=mean(job1),job2=mean(job2),job3=mean(job3),
                                                           job4=mean(job4),job5=mean(job5),job6=mean(job6),
                                                           job7=mean(job7))
          
          if(i==2) agg <- group_by(xc, par) %>% summarise( count=sum(count)/1000, 
                                                           wage=round(median(salary),0), age=median(age),
                                                           college=mean(college), female=mean(gender),
                                                           newbie=mean(year_hired),
                                                           job1=mean(job1),job2=mean(job2),job3=mean(job3),
                                                           job4=mean(job4),job5=mean(job5),job6=mean(job6),
                                                           job7=mean(job7))
          
          
          
          tab <- data.frame(t(as.matrix(agg)))
          tab[1,] <- round(tab[1,],0)
          tab[7,] <- round(tab[7,],0)
          tab <- cbind(names(agg),tab)
          tab <- tab[2:nrow(tab),2:5]
          
          vvv  <- sum(tab[1,])
          
          vec <- tab[1,]*100/vvv
          
          tab <- rbind(tab[1,],vec,tab[2:nrow(tab),])
          
          tab0 <- tab[1:2,]
          xxx <- 4
          tab1 <- tab[c(3,4,7),]
          tab2 <- tab[c(5,6,8:14),]
 
          tab0 <- format(round((tab0),1),nsmall=1) 
          tab2 <- format(round((tab2),2),nsmall=2) ; tab2
          tab1 <- format(round((tab1),0),nsmall=0) ; tab1
          
         tab <- rbind(tab0,tab1,tab2)
          
          
         
          if(i==1){tabx <- tab}else{tabx <- cbind(tabx,tab)}
          
        }
        
        
        
        cols <- ncol(tabx)*2+1
        res2 <- data.frame(matrix(0,nrow(tabx),cols))
    
        res2[,] <- "&"
        fill <- seq(from=2,to=cols,by=2)
        res2[,fill] <- tabx
        labeles <- c("Number ('000)","Share (pct)","Wages (median)",
                     "Age (median)","Start year","Bachelor's degree","Female",
                     "Low-skill Clerical","Social Worker","High Management","Middle Management","Low-skill Health",
                     "Teacher","Other")
        res2 <- cbind(labeles, res2)
        res2[,ncol(res2)] <- "\tabularnewline" ; res2[nrow(res2),ncol(res2)] <- "\tabularnewline[0.2cm]"
        print(res2, row.names=F)
     
        rm(summa1,summa2,res2,agg,ocs,xc,bureau,tab,tab0,tab1,tab2,tabx,vec)

      
    # =============================================================================================================================
          
      
    # Table C2 | Summary Statistics of Interviews  ================================================================================
        
        xc <- filter(cadunico,  !is.na(ocup) & salary>0 & education>=0   )
        
        xc$par <- ifelse(!is.na(xc$old_party),1,0)
        xc$par[!is.na(xc$new_party) & !is.na(xc$old_party)] <- 2
        xc$par[!is.na(xc$new_party) & is.na(xc$old_party)] <- 3
        xc$count <- 1
     
        agg <- group_by(xc, par) %>% summarise( count=sum(count)/1000, 
                                                enroll2009=mean(int_200912),  
                                                bf_elig=mean(int_elig), 
                                                bf_2012=mean(int_2012), 
                                                bf_new=mean(int_new),
                                                bf_exact1=mean(int_exact1),
                                                bf_exact2 =mean(int_exact2 )
                                                
        )
        agg$bf_elig_pc <- agg$bf_elig*100/agg$enroll2009
        agg$bf_2012_pc <- agg$bf_2012*100/agg$enroll2009
        agg$bbf_new_pc <- agg$bf_new*100/agg$enroll2009
        
        agg$bf_exact1_pc <- agg$bf_exact1*100/agg$enroll2009
        agg$bf_exact2_pc <- agg$bf_exact2*100/agg$enroll2009
        
        res <- data.frame(t(agg))
        res1 <- res[c(2,3,4,5,6,7,8),]
        res2 <- res[c(9:13),]
        
        res1 <- round(res1, 1)
        res2 <- round(res2, 1)
        
        resx <- rbind(res1,res2)
        rw <- c("Number of Interviewers ('000)","Total Interviews 2009-2012",
                "BF-eligible Interviews","Interviews in 2012","Interviews of New Entrants","Interviews at Income Threshold I","Interviews at Income Threshold II",
                "BF-eligible Interviews","Interviews in 2012","Interviews of New Entrants","Interviews at Income Threshold I","Interviews at Income Threshold II")
        a <- .LatexTable(resx,rw,spaces=c(2,4), dist=.2)
        print(a, row.names=FALSE)
        
        rm(a,agg,res,res1,res2,resx,xc) 
      
    # =============================================================================================================================
        
    
    # Table C3 | Balance ==========================================================================================================
 
        used <- MainData
        res <- data.frame(matrix(0,length(controls),5))
  
        
        for (i in 1:nrow(res)){
          
          used$var <- used %>% select(controls[i])
          
          bands <- rdbwselect(used$var, used$run, c = 0, p = 1, q = qq,bwselect = btype,   masspoints="off")$bws ; bands
          b1 <- bands[1]
          b2 <- bands[3]
           
          reg <- rdrobust(used$var,used$run,  b=b2, h=b1, q=qq, vce=vtype)
          reg2 <- rdrobust(used$var,used$run,  b=b2, h=b1, q=qq, level=90, vce=vtype )
          
          int <- round(reg$tau_cl[1],3)
          co <- round(reg$coef[1],3)
          se <- round(reg$se[3],3)
          ci1 <- round(reg2$ci[3,1],3)
          ci2 <- round(reg2$ci[3,2],3)
          n <- round(reg$N_h[1],0)
          stars <- ""
          if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
          if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
          
          res[i,1] <- paste(format(co,nsmall=3),stars,sep="")
          res[i,2] <- paste( "(" ,format(se,nsmall=3), ")",sep="")
          res[i,3] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
          
          res[i,4] <- format(int,nsmall=3)
          
          res[i,5] <- format(round(b1,2),nsmall=2)
          res[i,6] <- n
          
          print(i)
        }

     
        # Print for LaTex
        rw <- c("Mayor, age","Mayor, gender","Mayor, public servant","Mayor, college","Mayor, newcomer",
                "Mayor, PT","Mayor, PT's federal ally", "Mayor, PT's coalition",
                "Coalition size","Partisans, population","Partisans, bureaucracy","Partisans, interviewers")
        a <- .LatexTable(res,rw,spaces=c(1:11),dist=.2)
        print( a , row.names=FALSE)
      
        rm(used,a,res,reg,reg2,bands,b1,b2)
      
    # =============================================================================================================================
      
        
    # Table C4 | Robustness =======================================================================================================
    
        w <- MainData
        w$var1 <- w$outcome*100/w$total
        w$var2 <- w$outcome.2*100/w$total.2
        w$var3 <- w$outcome.3*100/w$total.3
        
        res <- data.frame(matrix(0,6,3))
   
        for (k in 1:4){
          pp <- 1
          qqq <- qq
          
          if(k == 1) {pp <- 2 ; qqq <- qq+1}
          region <- dummy_cols(w$reg)[,2:ncol(dummy_cols(w$reg))]
          parties <- dummy_cols(w$party)[,2:ncol(dummy_cols(w$party))]
          
          covars <- w %>% select(all_of(controls))
          if(k> 2) covars <- cbind(covars, region)
          if(k> 3) covars <- cbind(covars, parties)
          
          for (i in 1:3){
            
            if(i == 1) {w$var <- w$var1}
            if(i == 2) {w$var <- w$var2}
            if(i == 3) {w$var <- w$var3}
           
            if(k == 1) covars <- cbind(covars, parties, region)
            
            bands <- rdbwselect(w$var, w$run, c = 0, p = 1, q = qq, bwselect = btype, masspoints="off")$bws ; bands
            b1 <- bands[1]
            b2 <- bands[3]
            
            reg <- rdrobust(w$var,w$run,   b=b2, h=b1, q=qqq, p=pp , covs=covars) #; summary(reg)
            reg2 <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qqq, p=pp, covs=covars, level=90)
            reg3 <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qqq, p=pp )
            
            if(k > 1) reg <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq , covs=covars)
            if(k > 1) reg2 <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, level=90, covs=covars)
            
            int <- round(reg3$tau_cl[1],3)
            co <- round(reg$coef[1],3)
            se <- round(reg$se[3],3)
            ci1 <- round(reg$ci[3,1],3)
            ci2 <- round(reg$ci[3,2],3)
            n <- round(reg$N_h[1],0)
            stars <- ""
            if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
            if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
            
            res[1,i] <- paste(format(co,nsmall=3),stars,sep="")
            res[2,i] <- paste( "(" ,format(se,nsmall=3), ")",sep="")
            res[3,i] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
            
            res[4,i] <- format(int,nsmall=3)
            
            res[5,i] <- format(round(b1,2),nsmall=2)
            res[6,i] <- n
            
            
          }
          
          print(k)
          if(k==1){resx <- res}else{resx <- rbind(resx,res)}
          
        }
  
        # Print for LaTex
        rw <- c("Incumbency Effect"," ","C.I. (95pc)","Baseline")#,"Optimal Bandwidth","Municipalities")
        
        a <- .LatexTable(resx[7:10,],rw,spaces=c(4),dist=.2)
        print( a , row.names=FALSE)
        
        a <- .LatexTable(resx[13:16,],rw,spaces=c(4),dist=.2)
        print( a , row.names=FALSE)
        
        a <- .LatexTable(resx[19:22,],rw,spaces=c(4),dist=.2)
        print( a , row.names=FALSE)
    
        a <- .LatexTable(resx[1:4,],rw,spaces=c(4),dist=.2)
        print( a , row.names=FALSE)
        
        rw <- c("Bandwidth","Municipalities")
        
        a <- .LatexTable(resx[5:6,],rw,spaces=c(2),dist=.2)
        print( a , row.names=FALSE)
        
        rm(w,a,res,reg,reg2,bands,b1,b2,parties,region,covars,resx,reg3)
        
    # =============================================================================================================================
    
    
    # Table C5 | Mayor's Party ===================================================================================================
   
        partisans <- fread("data_partisans.csv")
        bureau <- fread("data_bureau.csv")
        
        wxx <- MainData
        wxx$var1 <- wxx$outcome*100/wxx$total
        wxx$var2 <- wxx$outcome.2*100/wxx$total.2
        wxx$var3 <- wxx$outcome.3*100/wxx$total.3
        
        newzy <- cadunico
        newzy$t <- newzy$coalition_side
        newzy$inter_rel <- ifelse(!is.na(newzy$new_party) & newzy$date_int < newzy$date_new,1,0)
        newzy <- filter(newzy, inter_rel==1)# &
        newzy <- filter(newzy,                      (newzy$new_party == newzy$party_win & newzy$t==1) |
                          (newzy$new_party == newzy$party_los & newzy$t==0)       )

        
        aa <- data.table(group_by(newzy,ibge,t) %>% summarise(outcome.new=sum(inter_rel)))
        wxx <- left_join(wxx,aa,by=c("ibge","t"))
        wxx[is.na(wxx)] <- 0
        
        newzy <- bureau
        newzy$t <- newzy$coalition_side
        newzy$inter_rel <- ifelse(!is.na(newzy$new_party) ,1,0)
        newzy <- filter(newzy, inter_rel==1 &
                          (
                            (newzy$new_party == newzy$party_win & newzy$t==1) |
                              (newzy$new_party == newzy$party_los & newzy$t==0)
                          ))

        aa <- data.table(group_by(newzy,ibge,t) %>% summarise(outcome2.new=sum(inter_rel)))
        wxx <- left_join(wxx,aa,by=c("ibge","t"))
        wxx[is.na(wxx)] <- 0
        
        newzy <- filter(partisans, late_term>=0 )
        newzy$t <- newzy$coalition_side
        aa <- data.table(group_by(newzy,ibge,t) %>% summarise(outcome3.new=sum(mayor)))
        wxx <- left_join(wxx,aa,by=c("ibge","t"))
        wxx[is.na(wxx)] <- 0
        wxx$outcome3.new <- wxx$outcome3.new - wxx$outcome2.new - wxx$outcome.new
        
        w <- wxx
        
        res <- data.frame(matrix(0,11,3))
        
        for (i in 1:3){
          
          if(i == 1) {w$var <- w$var1}
          if(i == 2) {w$var <- w$var2}
          if(i == 3) {w$var <- w$var3}
          
          bands <- rdbwselect(w$var, w$run, c = 0, p = 1, q = qq,bwselect = btype, masspoints="off")$bws ; bands
          b1 <- bands[1]
          b2 <- bands[3]
          
          reg <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, vce=vtype)
          reg2 <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, level=90 , vce=vtype)
           
          int <- round(reg$tau_cl[1],3)
          co <- round(reg$coef[1],3)
          se <- round(reg$se[3],3)
          ci1 <- round(reg$ci[3,1],3)
          ci2 <- round(reg$ci[3,2],3)
          n <- round(reg$N_h[1],0)
          stars <- ""
          if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
          if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
          
          col <- i
          
          res[1,col] <- paste(format(co,nsmall=3),stars,sep="")
          res[2,col] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
          
          res[3,col] <- format(int,nsmall=3)
          
          if(i == 1)  w$var <- w$outcome.new*100/w$total
          if(i == 2)  w$var <- w$outcome2.new*100/w$total.2
          if(i == 3)  w$var <- w$outcome3.new*100/w$total.3
          
          reg <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, vce=vtype)
          reg2 <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, level=90 , vce=vtype)
          ci1 <- round(reg$ci[3,1],3)
          ci2 <- round(reg$ci[3,2],3)
          
          stars <- ""
          if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
          if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
          
          co2 <- round(reg$coef[1],3)
          ci1 <- round(reg$ci[3,1],3)
          ci2 <- round(reg$ci[3,2],3)
          res[4,col] <- paste(format(co2,nsmall=3),stars,sep="")
          res[5,col] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
          res[6,col] <- paste(format(round(co2*100/co,0),nsmall=0),"%",sep="")
          
          if(i == 1)  w$var <-w$var1- w$outcome.new*100/w$total
          if(i == 2)  w$var <- w$var2- w$outcome2.new*100/w$total.2
          if(i == 3)  w$var <- w$var3- w$outcome3.new*100/w$total.3
          
          reg <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, vce=vtype)
          reg2 <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, level=90 , vce=vtype)
          ci1 <- round(reg$ci[3,1],3)
          ci2 <- round(reg$ci[3,2],3)
          
          stars <- ""
          if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
          if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
          
          co2 <- round(reg$coef[1],3)
          ci1 <- round(reg$ci[3,1],3)
          ci2 <- round(reg$ci[3,2],3)
          res[7,col] <- paste(format(co2,nsmall=3),stars,sep="")
          res[8,col] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
          res[9,col] <- paste(format(round(co2*100/co,0),nsmall=0),"%",sep="")
          
          
          res[10,col] <- format(round(b1,2),nsmall=2)
          res[11,col] <- n
          print(i)
          
        }
      
        rw <- c("Incumbency Effect","C.I. (95%)","Baseline",
                "Incumbency Effect (Outcome A)","C.I. (95%)","% of Total RD Effect",
                "Incumbency Effect (Outcome B)","C.I. (95%)","% of Total RD Effect",
                "Optimal Bandwidth","Number of Municipalities")
        a <- .LatexTable(res,rw,spaces=c(3,6,9,length(rw)), dist=.2)
        print(a, row.names=FALSE)
        
        rm(partisans,a,res,newzy,wxx,aa,w,reg,reg2,bands,bureau)
      
    # =============================================================================================================================
   
 
    # Table C6 | Electoral cycle ==================================================================================================
 
        partisans <- fread("data_partisans.csv")
        bureau <- fread("data_bureau.csv")
        
       
        wxx <- MainData
        wxx$var1 <- wxx$outcome*100/wxx$total
        wxx$var2 <- wxx$outcome.2*100/wxx$total.2
        wxx$var3 <- wxx$outcome.3*100/wxx$total.3
        
        newzy <- cadunico
        newzy$t <- newzy$coalition_side
        newzy$inter_rel <- ifelse(!is.na(newzy$new_party) & newzy$date_int < newzy$date_new,1,0)
        newzy <- filter(newzy, date_new < "2011-01-01"  )
         
        aa <- data.table(group_by(newzy,ibge,t) %>% summarise(outcome.new=sum(inter_rel)))
        wxx <- left_join(wxx,aa,by=c("ibge","t"))
        wxx[is.na(wxx)] <- 0
        
        newzy <- bureau
        newzy$t <- newzy$coalition_side
        newzy$inter_rel <- ifelse(!is.na(newzy$new_party) ,1,0)
        newzy <- filter(newzy, date_new < "2011-01-01"  )
         
        aa <- data.table(group_by(newzy,ibge,t) %>% summarise(outcome2.new=sum(inter_rel)))
        wxx <- left_join(wxx,aa,by=c("ibge","t"))
        wxx[is.na(wxx)] <- 0

        newzy <- filter(partisans, late_term==0 )
        newzy$t <- newzy$coalition_side
        aa <- data.table(group_by(newzy,ibge,t) %>% summarise(outcome3.new=sum(all)))
        wxx <- left_join(wxx,aa,by=c("ibge","t"))
        wxx[is.na(wxx)] <- 0
        wxx$outcome3.new <- wxx$outcome3.new - wxx$outcome2.new - wxx$outcome.new
        
        w <- wxx
         
        res <- data.frame(matrix(0,11,3))
        
        for (i in 1:3){
          
          if(i == 1) {w$var <- w$var1}
          if(i == 2) {w$var <- w$var2}
          if(i == 3) {w$var <- w$var3}
          
          bands <- rdbwselect(w$var, w$run, c = 0, p = 1, q = qq,bwselect = btype, masspoints="off")$bws ; bands
          b1 <- bands[1]
          b2 <- bands[3]
          
          reg <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, vce=vtype)
          reg2 <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, level=90 , vce=vtype)
         
          int <- round(reg$tau_cl[1],3)
          co <- round(reg$coef[1],3)
          se <- round(reg$se[3],3)
          ci1 <- round(reg$ci[3,1],3)
          ci2 <- round(reg$ci[3,2],3)
          n <- round(reg$N_h[1],0)
          stars <- ""
          if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
          if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
          
          col <- i
          
          res[1,col] <- paste(format(co,nsmall=3),stars,sep="")
          res[2,col] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
          
          res[3,col] <- format(int,nsmall=3)
          
          if(i == 1)  w$var <- w$outcome.new*100/w$total
          if(i == 2)  w$var <- w$outcome2.new*100/w$total.2
          if(i == 3)  w$var <- w$outcome3.new*100/w$total.3
          
          reg <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, vce=vtype)
          reg2 <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, level=90 , vce=vtype)
          ci1 <- round(reg$ci[3,1],3)
          ci2 <- round(reg$ci[3,2],3)
     
          stars <- ""
          if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
          if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
          
          co2 <- round(reg$coef[1],3)
          ci1 <- round(reg$ci[3,1],3)
          ci2 <- round(reg$ci[3,2],3)
          res[4,col] <- paste(format(co2,nsmall=3),stars,sep="")
          res[5,col] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
          res[6,col] <- paste(format(round(co2*100/co,0),nsmall=0),"%",sep="")
          
          if(i == 1)  w$var <-w$var1- w$outcome.new*100/w$total
          if(i == 2)  w$var <- w$var2- w$outcome2.new*100/w$total.2
          if(i == 3)  w$var <- w$var3- w$outcome3.new*100/w$total.3
          
          reg <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, vce=vtype)
          reg2 <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, level=90 , vce=vtype)
          ci1 <- round(reg$ci[3,1],3)
          ci2 <- round(reg$ci[3,2],3)
       
          stars <- ""
          if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
          if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
          
          co2 <- round(reg$coef[1],3)
          ci1 <- round(reg$ci[3,1],3)
          ci2 <- round(reg$ci[3,2],3)
          res[7,col] <- paste(format(co2,nsmall=3),stars,sep="")
          res[8,col] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
          res[9,col] <- paste(format(round(co2*100/co,0),nsmall=0),"%",sep="")
          
          
          res[10,col] <- format(round(b1,2),nsmall=2)
          res[11,col] <- n
          print(i)
          
        }

   
        rw <- c("Incumbency Effect","C.I. (95%)","Baseline",
                "Incumbency Effect (Outcome A)","C.I. (95%)","% of Total RD Effect",
                "Incumbency Effect (Outcome B)","C.I. (95%)","% of Total RD Effect",
                "Optimal Bandwidth","Number of Municipalities")
        a <- .LatexTable(res,rw,spaces=c(3,6,9,length(rw)), dist=.2)
        print(a, row.names=FALSE)
        
        rm(partisans,a,res,newzy,wxx,aa,w,reg,reg2,bands,bureau)
        
      
    # =============================================================================================================================

  
    # Table C7 | Comparison | CU Expansion ========================================================================================
      
        
        w <- MainData
       
        # Create Table
        res <- data.frame(matrix(0,8,3))
        for (i in 1:3){
          
          if(i == 1) {w$var <- w$outcome*100/w$total}
          if(i == 2) {w$var <- w$outcome.2*100/w$total.2}
          if(i == 3) {w$var <- w$outcome.3*100/w$total.3}
          w$i <-    (w$pot_enrollment-mean(w$pot_enrollment))/sd(w$pot_enrollment)
          
          # Set bands
          bands <- rdbwselect(w$var, w$run, c = 0, p = 1, q = qq,bwselect = btype, masspoints="off")$bws ; bands
          b1 <- bands[1]
          b2 <- bands[3]
          
          data <- filter(w,  abs(run) < b1   )
          data$wght <- b1 - abs(data$run)
          regols <- felm( var ~  t*run *i    +age+gender+public_servant+newcomer+education
                          +past_memb1+past_memb2+past_memb3+pt_coalition+size
                          | party+reg | 0 | ibge ,weights=data$wght , data=data) 
          rr <- summary(regols)$coefficients 
          vv <- vcov(regols)
       
          co <- round(rr["t",1],3)
          se <- round(rr["t",2],3)
         
          stars <- ""
          if(rr["t",4]<=0.10) stars <- "+"
          if(rr["t",4]<=0.05) stars <- "*"
          
          res[1,i] <- paste(format(co,nsmall=3),stars,sep="")
          res[2,i] <- paste( "(" ,format(se,nsmall=3), ")",sep="")
          
          co <- round(rr["i",1],3)
          se <- round(rr["i",2],3)
          
          stars <- ""
          if(rr["i",4]<=0.10) stars <- "+"
          if(rr["i",4]<=0.05) stars <- "*"
          
          res[3,i] <- paste(format(co,nsmall=3),stars,sep="")
          res[4,i] <- paste( "(" ,format(se,nsmall=3), ")",sep="")
          
          
          co <- round(rr["t:i",1],3)
          se <- round(rr["t:i",2],3)
          
          stars <- ""
          if(rr["t:i",4]<=0.10) stars <- "+"
          if(rr["t:i",4]<=0.05) stars <- "*"
          
          res[5,i] <- paste(format(co,nsmall=3),stars,sep="")
          res[6,i] <- paste( "(" ,format(se,nsmall=3), ")",sep="")

          
          res[7,i] <- format(round(b1,2),nsmall=2)
          res[8,i] <- length(unique(data$ibge))
          print(i)
          
        }
        
        # Print Table in LaTex format
        rw <- c("Incumbency Effect (A)"," ","Potential Enrollment (B)"," ","(A) x (B)", "","Optimal Bandwidth",
                "Municipalities")
        a <- .LatexTable(res,rw,spaces=c(2,4,6,8),dist=.2)
        print( a , row.names=FALSE)
        
        rm(w,a,res,regols,rr,vv,co,se,stars,data,bands)
        
        
    # =============================================================================================================================
     
  
    # Table C8 | Balance Matched ==================================================================================================
      
      # get data on all individuals   
      rai <- fread("data_bureau.csv")
      cad <- cadunico
  
      # get data on matched individuals
      m <- fread("data_matched.csv")
      m0 <- m[,c(3:11,2)]
      m1 <- m[,c(12:20,2)]
      
      names(m0) <- names(m1) <- c("admin","social", "health", "lowskill", "manag", "salary", "education", "year_born", "year_hired", "gender")
      
      cad <- filter(cad,  !is.na(ocup) & salary>0 & education>=0   )
      rai <- filter(rai, !is.na(ocup) & salary>0 & education>=0   )
      
      # Create variables
      rai$other <- ifelse(rai$admin+rai$social+rai$health+rai$lowskill+rai$manag==0,1,0)
      cad$other <- ifelse(cad$admin+cad$social+cad$health+cad$lowskill+cad$manag==0,1,0)
      m0$other <- ifelse(m0$admin+m0$social+m0$health+m0$lowskill+m0$manag==0,1,0)
      m1$other <- ifelse(m1$admin+m1$social+m1$health+m1$lowskill+m1$manag==0,1,0)
      
      rai$col <- ifelse(rai$education > 8,1,0)
      cad$col <- ifelse(cad$education > 8,1,0)
      m1$col <- ifelse(m1$education > 8,1,0)
      m0$col <- ifelse(m0$education > 8,1,0)
      
      rai$age <- 2008-rai$year_born
      cad$age <- 2008-cad$year_born
      m0$age <- 2008-m0$year_born
      m1$age <- 2008-m1$year_born
      
      rai$yin2 <- 2008-rai$year_hired
      cad$yin2 <- 2008-cad$year_hired
      m0$yin2 <- 2008-m0$year_hired
      m1$yin2 <- 2008-m1$year_hired
      

      var_table <- c("admin","social","health","lowskill","manag","other","gender","age","education","yin2","salary")
 
      #Run averages
      matched <- data.frame(matrix(0,length(var_table),4))
      for (i in 1:length(var_table)){
        
        roundi <- 3
        if(i > 10) roundi <- 1
        
        cad$v <- cad %>% select(var_table[i])
        rai$v <- rai %>% select(var_table[i])
        m0$v <- m0 %>% select(var_table[i])
        m1$v <- m1 %>% select(var_table[i])
        
        tt <- t.test(cad$v ,rai$v, na.rm=TRUE )
        pv <- tt$p.value
        star <- ""
        if(pv<.1) star <- "+"
        if(pv<.05) star <- "*"
        
        tt <- t.test(m0$v ,m1$v, na.rm=TRUE )
        pv <- tt$p.value
        star1 <- ""
        if(pv<.1) star1 <- "+"
        if(pv<.05) star1 <- "*"
        
        
        matched[i,1] <- format(round(mean(cad$v, na.rm=TRUE),roundi),nsmall=roundi)
        matched[i,2] <- paste0(format(round(mean(rai$v, na.rm=TRUE)-mean(cad$v , na.rm=TRUE),roundi),nsmall=roundi),star)
        
        matched[i,3] <- format(round(mean(m0$v, na.rm=TRUE),roundi),nsmall=roundi)
        
        if(i <= 7){        matched[i,4] <- "-" }else{
        matched[i,4] <- paste0(format(round(mean(m1$v, na.rm=TRUE)-mean(m0$v , na.rm=TRUE),roundi),nsmall=roundi),star1)
        }
      }   
      mmm <- c(nrow(cad),nrow(rai),nrow(m),"-")
      matched <- rbind(matched,mmm)
      
      # Print for LaTex
      rw <- c("Low-skill Clerical","High-skill SS", "Health Workers and Teachers","Low-skill General", 
              "High Management","Other Occupation", "Gender","Age","Education","Years Employed",  "Wage","Observations")
      a <- .LatexTable(matched,rw,spaces=c(11,12), dist=.2)
      print( a , row.names=FALSE)
      
      rm(a,cad,rai,m0,m1,mmm,m,matched,tt)

    # =============================================================================================================================
 
      
    # Table C9 | teachers and health workers ======================================================================================
      
      w1 <- w2 <- MainData
      w1$i <- 0
      w2$i <- 1
      wxx <- rbind(w1,w2)
      
      newzy <- fread("data_bureau.csv")
      newzy$t <- newzy$coalition_side
      newzy$ocup2 <- substring(newzy$ocup,1,3)
      newzy$inter_tot <- 1 
      newzy$inter_rel <- ifelse(!is.na(newzy$new_party) ,1,0)
      newzy$i <- ifelse(newzy$ocup2 %in% c(515) | newzy$ocup %in% c(352210),1,0 ) #17% of bureaucrats 
      
      
      aa <- data.table(group_by(newzy,ibge,t,i) %>% summarise(outcome.new=sum(inter_rel),inter_tot=sum(inter_tot)))
      bb <- data.table(group_by(aa,ibge,i) %>% summarise(total.new=sum(inter_tot)))
      aa$inter_tot <- NULL
      
      wxx <- left_join(wxx,bb,by=c("ibge","i"))
      
      wxx <- left_join(wxx,aa,by=c("ibge","t","i"))
      wxx[is.na(wxx)] <- 0
      
      wxx$var <- wxx$outcome.new*100 / wxx$total.new
      wxx <- filter(wxx, !is.na(var))
      
      bands <- rdbwselect(wxx$var, wxx$run, c = 0, p = 1, q = qq,bwselect = btype , masspoints="off")$bws 
      b1 <- bands[1] 
      b2 <- bands[3]
      
      res <- data.frame(matrix(0,5,3))
      j <- 1
      for (j in 0:1){
        
        wx <- filter(wxx, i == j)  
         
        if(j==0) b1a <- b1
        if(j==1) b1b <- b1
        
        covars <- wx %>% select(all_of(controls))
        covars <- cbind(covars)
        
        reg <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars,vce=vtype) 
        reg2 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars, level=90, vce=vtype) 
        reg3 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, vce=vtype) 
        
        if(j==0) c1 <- reg$coef[1]
        if(j==1) c2 <- reg$coef[1]
        co <- round(reg$coef[1],3)
        se <- round(reg$se[3],3)
        ci1 <- round(reg$ci[3,1],3)
        ci2 <- round(reg$ci[3,2],3)
        n <- round(reg$N_h[1],0)
        stars <- ""
        if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
        if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
        bas <- round(reg3$tau_cl[1],3)
        if(j==0) bas1 <- bas
        if(j==1) bas2 <- bas
        
        res[1,j+1] <- paste(format(co,nsmall=3),stars,sep="")
        #res[2,i+1] <- paste( "(" ,format(se,nsmall=3), ")",sep="")
        res[2,j+1] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
        res[3,j+1] <- paste(format(bas,nsmall=3),sep="")
        
        dt <- filter(wx, abs(run) < b1 )
        res[4,j+1] <- length(unique(dt$ibge))
        res[5,j+1] <- format(round(b1,2),nsmall=2) 
      }
      
      
      cities <- wxx
      cities <- unique(cities, by="ibge")
      cities$x <- 1
      cities <- cities %>% select(ibge,x)
      
      set.seed(14627)
      results <- boot(data=cities, base=wxx, statistic=.DiffRDD1 , R=500,con=1,v1="var",v2="var",
                      b1a=b1a, b1b=b1b,  parallel = "snow", ncpus = 4, byvar=1)
      
      ci = boot.ci(results, type=c("perc","basic"), conf=.95, vce=vtype)  
      ci2 = boot.ci(results, type=c("perc","basic"), conf=.90, vce=vtype) 
      
      stars <- ""
      if(ci2$perc[4]*ci2$perc[5]>0) stars <- "+"
      if(ci$perc[4]*ci$perc[5]>0) stars <- "*"
      
      ci1 <- round(ci$perc[4],3)
      ci2 <- round(ci$perc[5],3)
      
      res[1,3] <- paste(format(round(c1-c2,3),nsmall=3),stars,sep="")
      res[2,3] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
      res[3,3] <- paste(format(round(bas1-bas2,3),nsmall=3),sep="")
      res[4:5,3] <- "-"
      res1 <- res
      ##
      
      
      rw <- c("Incumbency Effect","C.I. (95%)","Baseline","Municipalities","Optimal Bandwidth")
      a <- .LatexTable(res1,rw,spaces=c(5),dist=0.2)
      print(a, row.names=FALSE)
      
      rm(wxx,wx,w1,w2,reg,res,bands,b1,b2,covars,newzy,aa,bb,reg2,reg3,res1,results,ci,cities,dt,a)
      
      
    # =============================================================================================================================
      
  
    # Table C10 | demographic effects =============================================================================================
        
        w1 <- w2 <- MainData
        w1$i <- 0
        w2$i <- 1
        v <- rbind(w1,w2)
        v$var1 <- v$outcome*100/v$total
        
        # Get individual interviewers
        c <- cadunico 
        c$t <- c$coalition_side
        c <- data.table( c %>% group_by(ibge) %>% mutate(med_salary=median(salary)))
        
        # Set bands
        bands <- rdbwselect(v$var1, v$run, c = 0, p = 1, q = qq,bwselect = btype, masspoints="off")$bws ; bands
        b1 <- bands[1]
        b2 <- bands[3]
          
        # Create Table
        for (i in 1:6){
          
            c$inter_tot <- 1 
            c$inter_rel <- ifelse(!is.na(c$new_party) & as.numeric(c$date_int) < (as.numeric(c$date_new))  ,1,0)
            
            if(i==1) c$i <- c$civil 
            if(i==2) c$i <- ifelse(c$education >= 8,1,0 )
            if(i==3) c$i <- ifelse(c$gender == 1,1,0 )
            if(i==4) c$i <- ifelse(c$salary < c$med_salary,1,0 )
            if(i==5) c$i <- ifelse(c$year_hired >= 2005,1,0 )
            if(i==6) c$i <- ifelse(!is.na(c$old_party),1,0 )
            
            aa <- data.table(group_by(c,ibge,t,i) %>% summarise(outcome.new=sum(inter_rel),inter_tot=sum(inter_tot)))
            bb <- data.table(group_by(aa,ibge,i) %>% summarise(total.new=sum(inter_tot)))
            aa$inter_tot <- NULL
            
            newby <- left_join(v,aa,by=c("ibge","t","i"))
            newby[is.na(newby)] <- 0
            newby <- left_join(newby,bb,by=c("ibge","i"))
            newby[is.na(newby)] <- 0
            
            
            newby$var <- newby$outcome.new*100 / newby$total.new
            
            res <- data.frame(matrix(0,3,3))
            
            for (j in 0:1){
              
              
              wx <- filter(newby, i == j)  
               
              covars <- wx %>% select(age,gender,newcomer,size, public_servant,education,past_memb1,past_memb2,past_memb3,pt_coalition,pt_mayor,pt_federal)
              covars <- cbind(covars)
              
              reg <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars) 
              reg2 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars, level=90) 
              reg3 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq) 
              
              if(j==0) c1 <- reg$coef[1]
              if(j==1) c2 <- reg$coef[1]
              co <- round(reg$coef[1],3)
              se <- round(reg$se[3],3)
              ci1 <- round(reg$ci[3,1],3)
              ci2 <- round(reg$ci[3,2],3)
              
              bas1 <- round(reg3$tau_cl[1],3)
              if(j==0) bas1a <- bas1
              if(j==1) bas1b <- bas1
              
              
              
              n <- round(reg$N_h[1],0)
              stars <- ""
              if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
              if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
              
              res[1,j+1] <- paste(format(co,nsmall=3),stars,sep="")
              res[2,j+1] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
              dt <- filter(wx, abs(run) < b1 )
          
              res[3,j+1] <- format( round(reg$tau_cl[1],3) ,nsmall=3)
              ex1 <- length(unique(dt$ibge))
              ex2 <- format(round(b1,2),nsmall=2) 
            }
            
            
            
            wx <- newby
            cities <- unique(wx, by="ibge")
            cities <- cities %>% select(ibge)
            
            set.seed(14627)
            results <- boot(data=cities, base=wx, statistic=.DiffRDD1 , R=500,con=1,byvar=1,v1="var",v2="var",
                            b1a=b1a, b1b=b1b,
                            parallel = "snow", ncpus = 4)
            
            ci = boot.ci(results, type=c("perc","basic"), conf=.95)  ;  ci 
            ci2 = boot.ci(results, type=c("perc","basic"), conf=.90) ;  ci2
            
            stars <- ""
            if(ci2$perc[4]*ci2$perc[5]>0) stars <- "+"
            if(ci$perc[4]*ci$perc[5]>0) stars <- "*"
            
            ci1 <- round(ci$perc[4],3)
            ci2 <- round(ci$perc[5],3)
            
            res[1,3] <- paste(format(round(c1-c2,3),nsmall=3),stars,sep="")
            res[2,3] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
            res[3,3] <- format(round(bas1a-bas1b,3),nsmall=3)
            if(i==1) res1 <- res
            if(i==2) res2 <- res
            if(i==3) res3 <- res
            if(i==4) res4 <- res
            if(i==5) res5 <- res
            if(i==6) res6 <- res
            print(i)
            
          }
          
       
        # Print Table Panels in LaTex format  
          rw <- c("Incumbency Effect","C.I. (95%)"," Baseline")
          a <- .LatexTable(res1,rw,spaces=c(length(rw)),dist=.2)
          print(a, row.names=FALSE)
          
          a <- .LatexTable(res2,rw,spaces=c(length(rw)),dist=.2)
          print(a, row.names=FALSE)
          
          a <- .LatexTable(res3,rw,spaces=c(2,length(rw)),dist=.2)
          print(a, row.names=FALSE)
          
          a <- .LatexTable(res4,rw,spaces=c(2,length(rw)),dist=.2)
          print(a, row.names=FALSE)
          
          a <- .LatexTable(res5,rw,spaces=c(2,length(rw)),dist=.2)
          print(a, row.names=FALSE)
          
          a <- .LatexTable(res6,rw,spaces=c(2,length(rw)),dist=.2)
          print(a, row.names=FALSE)
          
          rm(v,wx,w1,w2,reg,res,bands,b1,b2,covars,newby,aa,bb,reg2,reg3,res1,res2,res3,res4,res5,res6,results,ci,cities,dt,a,c)
          gc()
          
          
    # =============================================================================================================================
        
 
    # Table C11 | Balance Interviewers Incumbent vs Opposition ====================================================================
      
      # Create Variables
      xxx <- filter(cadunico, coalition_side>=0 & !is.na(remains_bureaucrat) )
      xxx$t <- xxx$coalition_side
      xxx$v1 <- ((xxx$salary)/100)
      xxx$boss1 <- ifelse(xxx$ocup < 192000 ,1,0)
      xxx$i <- xxx$civil

      # Run Averages
      res <- data.frame(matrix(0,8,3))
      for (i in 1:8){
        
        if(i==1) xxx$var <- xxx$v1
        if(i==2) xxx$var <- xxx$i
        if(i==3) xxx$var <- xxx$boss1
        if(i==4) xxx$var <- 2008-xxx$year_hired
        if(i==5) xxx$var <- 2008-xxx$year_born
        if(i==6) xxx$var <- xxx$gender
        if(i==7) xxx$var <- xxx$education
        if(i==8) xxx$var <- xxx$education12
  
        reg <- felm(var ~ t        |  0 | 0 | ibge , data=xxx) ; r <- summary(reg)$coefficients
        
        res[i,1] <- r["t",1]
        res[i,2] <- r["t",2]
        res[i,3] <- r["(Intercept)",1]
        
      }
      
      res$st <- ""
      res$st[abs(res$X1/res$X2)>=1.65] <- "+"
      res$st[abs(res$X1/res$X2)>=1.96] <- "*"
      
      res$baseline <- paste0(format(round(res$X3,3),nmall=3))
      res$effect <- paste0(format(round(res$X1,3),nmall=3),res$st)
      res$sd <- paste0("(",format(round(res$X2,3),nmall=3),")")
      
      # Print for LaTex
      res2 <- res %>% select(baseline,effect,sd)
      rw <- c("Salary",
              "Job Security",
              "Managerial Position",
              "Seniority",
              "Age",
              "Gender",
              "Education",
              "Education (2012)+")
      a <- .LatexTable(res2,rw,spaces=c(length(rw)), dist=0.2)
      print(a, row.names=FALSE)
 
      rm(a,r,reg,res,res2,xxx)
      
    # =============================================================================================================================
      
      
    # Table C12 | 2nd vs 3rd Place ================================================================================================
        
        w <- filter(MainData, !is.na(mv3rd))
        w$run_old <- w$run
        w$run <- (w$t*2-1)*w$mv3rd
 
        res <- data.frame(matrix(0,6,3))
        w$var1 <- w$outcome.3rd*100/w$total
        w$var2 <- w$outcome2.3rd*100/w$total.2
        w$var3 <- (w$outcome3.3rd-w$outcome2.3rd-w$outcome.3rd)*100/w$total.3

        for (i in 1:3){
          
          if(i == 1) {w$var <- w$var1}
          if(i == 2) {w$var <- w$var2}
          if(i == 3) {w$var <- w$var3}
          
          bands <- rdbwselect(w$var, w$run, c = 0, p = 1, q = qq,bwselect = btype, masspoints="off")$bws 
          b1 <- bands[1]
          b2 <- bands[3]

          
          covars <- w %>% select(age,gender,newcomer, size,public_servant,education,past_memb1,past_memb2,past_memb3,pt_coalition,pt_mayor,pt_federal)
          covars <- cbind(covars)
          
          reg <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, vce=vtype, covs=covars)
          reg2 <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, level=90 , vce=vtype , covs=covars)
          reg3 <- rdrobust(w$var,w$run,  b=b2, h=b1, q=qq, vce=vtype)
          
          
          int <- round(reg3$tau_cl[1],3)
          co <- round(reg$coef[1],3)
          se <- round(reg$se[3],3)
          ci1 <- round(reg$ci[3,1],3)
          ci2 <- round(reg$ci[3,2],3)
          n <- round(reg$N_h[1],0)
          stars <- ""
          if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
          if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
          
          res[1,i] <- paste(format(co,nsmall=3),stars,sep="")
          res[2,i] <- paste( "(" ,format(se,nsmall=3), ")",sep="")
          res[3,i] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
          
          res[4,i] <- format(int,nsmall=3)
          
          res[5,i] <- format(round(b1,2),nsmall=2)
          res[6,i] <- n
          print(i)
          
        }
        
        rw <- c("Runner-up Effect"," ","C.I. (95%)","Baseline","Optimal Bandwidth","Municipalities")
        a <- .LatexTable(res,rw,spaces=c(5), dist=.2)
        print(a, row.names=FALSE)
        
        rm(a,bands,covars,reg,reg2,reg3,res,w)
        
    
        
    # =============================================================================================================================
          
          
    # Table E1 | Competition and Reelection =======================================================================================
          
          c <- MainData
          cc <- filter(c, inc==1 & t == 1)
       
          for (i in 1:3){
            
            c$var <- c$outcome*100/c$total
            bands <- rdbwselect(c$var, c$run, c = 0, p = 1, q = qq,bwselect = btype, masspoints="off")$bws ; bands
            b1 <- bands[1]
            b2 <- bands[3]
            
            if(i==1) c$i <- ifelse(c$hh < median(c$hh[abs(c$run)<b1]),1,0 )
            if(i==2) c$i <- ifelse(c$cand_tot >= 3,1,0 )
            if(i==3) c$i <- ifelse(c$ibge %in% cc$ibge,0,1 )
            
            
            res <- data.frame(matrix(0,5,3))
            
            for (j in 0:1){
              
              wx <- filter(c, i == j)  
              covars <- wx %>% select(age,gender,newcomer, size,public_servant,education,past_memb1,past_memb2,past_memb3,pt_coalition,pt_mayor,pt_federal)
              covars <- cbind(covars)
              
              reg <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, covs=covars) 
              reg2 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq, level=90, covs=covars) 
              reg3 <- rdrobust(wx$var,wx$run,b=b2, h=b1, q=qq) 
              
              if(j==0) c1 <- reg$coef[1]
              if(j==1) c2 <- reg$coef[1]
              co <- round(reg$coef[1],3)
              se <- round(reg$se[3],3)
              ci1 <- round(reg$ci[3,1],3)
              ci2 <- round(reg$ci[3,2],3)
              
              bas1 <- round(reg3$tau_cl[1],3)
              if(j==0) bas1a <- bas1
              if(j==1) bas1b <- bas1
              
              
              
              n <- round(reg$N_h[1],0)
              stars <- ""
              if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
              if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
              
              res[1,j+1] <- paste(format(co,nsmall=3),stars,sep="")
              res[2,j+1] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
              dt <- filter(wx, abs(run) < b1 )
              
              res[3,j+1] <- format( round(bas1,3) ,nsmall=3)
              res[4,j+1] <- length(unique(dt$ibge))
              res[5,j+1] <- format(round(b1,2),nsmall=2) 
            }
            
            res
            
            wx <- c
            cities <- unique(wx, by="ibge")
            cities <- cities %>% select(ibge)
            set.seed(14627)
            results <- boot(data=cities, base=wx, statistic=.DiffRDD1 , R=500,con=1,byvar=1,v1="var",v2="var",
                            b1a=b1, b1b=b1,parallel = "snow", ncpus = 4)
            
            ci = boot.ci(results, type=c("perc","basic"), conf=.95)  ;  ci 
            ci2 = boot.ci(results, type=c("perc","basic"), conf=.90) ;  ci2
            
            stars <- ""
            if(ci2$perc[4]*ci2$perc[5]>0) stars <- "+"
            if(ci$perc[4]*ci$perc[5]>0) stars <- "*"
            
            ci1 <- round(ci$perc[4],3)
            ci2 <- round(ci$perc[5],3)
            
            res[1,3] <- paste(format(round(c1-c2,3),nsmall=3),stars,sep="")
            res[2,3] <- paste( "[" ,format(ci1,nsmall=3),",",format(ci2,nsmall=3), "]",sep="")
            res[3,3] <- format(round(bas1a-bas1b,3),nsmall=3)
            res[4,3] <- format(round(as.numeric(res[4,1])-as.numeric(res[4,2]),0),nsmall=0)
            res[5,3] <- " "
            
            if(i==1) res1 <- res
            if(i==2) res2 <- res
            if(i==3) res3 <- res
            if(i==4) res4 <- res
      
            print(i)
            
          }
          
     
          rw <- c("Incumbency Effect","C.I. (95%)"," Baseline","Municipalities","Bandwidth")
          a <- .LatexTable(res1,rw,spaces=c(length(rw)),dist=.2)
          print(a, row.names=FALSE)
          
          a <- .LatexTable(res2,rw,spaces=c(length(rw)),dist=.2)
          print(a, row.names=FALSE)
          
          a <- .LatexTable(res3,rw,spaces=c(length(rw)),dist=.2)
          print(a, row.names=FALSE)
          
          rm(a,c,bands,cc,res,wx,res1,res2,res3,ci,cities,covars,dt,reg,reg2,reg3,results)

    # =============================================================================================================================
          
         
          
    ## F I G U R E S           
          
    
    # Figure C1 | RDDs =============================================================================================================
    
      w <- MainData
      w$var1 <- w$outcome*100/w$total
      w$var2 <- w$outcome.2*100/w$total.2
      w$var3 <- w$outcome.3*100/w$total.3
          
      band <- 15
      binss <- 25
      confi <- .95
      scl1 <- 7
      scl2 <- 3
      data <- filter(w, abs(run) < band   )
      
      data$var <- data$var1
      sequ <- seq(from=-band, to=band, length=7)
      limi1 <- 0 ; limi2 <- 6
      rd1 <- .RDPlot(data=data,cut=0,band=band,bins=binss,degree=1,"","% of Interviewers","Interviewers",
                     con=confi,sz=.5,brk=sequ, limi1,limi2) ; rd1
      
      
      data$var <- data$var2
      sequ <- seq(from=-band, to=band, length=7)
      limi1 <- 0 ; limi2 <- 3
      rd2 <- .RDPlot(data=data,cut=0,band=band,bins=25,degree=1,"","% of Other Bureaucrats","Other Bureaucrats",
                     con=confi,sz=.5,brk=sequ, limi1,limi2) ; rd2
      
      
      data$var <- data$var3 
      sequ <- seq(from=-band, to=band, length=7)
      limi1 <- 0 ; limi2 <- 1.5
      rd3 <- .RDPlot(data=data,cut=0,band=band,bins=binss,degree=1,"","% of Voters","Voters",
                     con=confi,sz=.5,brk=sequ, limi1,limi2) ; rd3
      
      grid.arrange(rd1, rd2, rd3, nrow = 3, bottom="Margin of Victory (%)")
      
      rm(w,band,binss,confi,scl1,scl2,data,sequ,limi1,limi2,rd1,rd2,rd3)
    
    # =============================================================================================================================
   

    # Figure C2 | Bandwidth Plots ====================================================================================================
    
      w <- MainData
      w$var1 <- w$outcome*100/w$total
      w$var2 <- w$outcome.2*100/w$total.2
      w$var3 <- w$outcome.3*100/w$total.3
      
      n <- 10
      nn <- (1:n)/5
   
      for (qw in 1:3){
        
        mat <- data.table(matrix(0,n,7))
        if(qw == 1) {w$var <- w$var1}
        if(qw == 2) {w$var <- w$var2}
        if(qw == 3) {w$var <- w$var3}
        
        if(qw == 1) {l1 <- -2 ; l2 <- 6}
        if(qw == 2) {l1 <- -.4 ; l2 <- .8}
        if(qw == 3) {l1 <- -.4 ; l2 <- .4}
        
        
        bands <- rdbwselect(w$var, w$run, c = 0, p = 1, q = qq, bwselect = btype, masspoints="off")$bws 
        b1 <- bands[1]
        b2 <- bands[3]
        
        for (i in 1:n){
          a <- nn[i]
          reg <- rdrobust(w$var, w$run, p=1,q=qq ,h=b1*a, b=b2*a ) 
          reg2 <- rdrobust(w$var, w$run, p=1,q=qq ,h=b1*a, b=b2*a, level=90 ) 
          mat[i,1] <- a*b1
          mat[i,3] <- reg$ci[3,1]
          mat[i,4] <- reg$ci[3,2]
          mat[i,5] <- reg2$ci[3,1]
          mat[i,6] <- reg2$ci[3,2]
          
          mat[i,2] <- reg$coef[1]
          mat[i,7] <- a*100
        }
        mat
        
        if(qw == 1) p.title <- "CadUnico Interviewers"
        if(qw == 2) p.title <- "Other Bureaucrats"
        if(qw == 3) p.title <- "Voters"
        
        names(mat) <- c("band","co","ci95a","ci95b","ci90a","ci90b","band2")
        pl <- ggplot(data=mat, aes(x=band2, y=co))   +    theme_bw() +
          ylim(l1,l2)+
          theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
          geom_point(aes(x=band2, y=co), colour="gray15",  size=4, alpha=.75) +
          geom_hline(yintercept = 0 , lty=3, alpha=.75) +
          theme(plot.title = element_text(family = "A", size=11))+
          xlab("% of Optimal Bandwidth") + ylab("RDD Effect") +
          theme(axis.text = element_text(colour="black", size=11, family="A")) +
          theme(axis.title = element_blank()) +
          theme(plot.title = element_text(colour="black", size=11, family="A")) +
          geom_errorbar(aes(x=band2, ymin=ci95a, ymax=ci95b), colour="gray15",  size=1, alpha=.75, width=0) +
          geom_errorbar(aes(x=band2, ymin=ci90a, ymax=ci90b), colour="gray25",  size=3, alpha=.5, width=0) +
          scale_x_continuous(breaks=nn*100)+
          ggtitle(p.title)
        
        if(qw == 1) pl1 <- pl
        if(qw == 2) pl2 <- pl
        if(qw == 3) pl3 <- pl
        print(qw)
        
      }
      grid.arrange(pl1, pl2, pl3, nrow = 3, bottom="Bandwidth (percentage of optimal)")
      
      rm(w,n,nn,qw,pl,pl1,pl2,pl3,bands,b1,b2,mat,reg,reg2)
    
    # =============================================================================================================================
    
          
    # Figure C3 | Histogram Potential enrollment ==================================================================================
        
        used <- filter(MainData,t==1)
       
        ggplot(used, aes(x=pot_enrollment)) + theme_bw() +
          theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
          geom_histogram(bins=50,fill="gray75",color="gray15") +
          theme(axis.title = element_text(colour="black", size=12, family="A")) +
          theme(axis.title.y = element_blank()) +xlab("Potential for CadUnico Enrollments (share of population)")+
          theme(axis.text = element_text(colour="black", size=12, family="A")) 
        
        rm(used)
    
    # =============================================================================================================================
        
              
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
       